//////////////////////////////////////////////////////////////////////////////////
//  PVR Texture Compressor														//
//  Author: Chirantan Ekbote													//
//	Mailto: ekbote.1@osu.edu													//
//																				//
//  Harvard School of Engineering and Applied Sciences							//
//  The Ohio State University													//
//																				//
//	PVRTexCUDA.cu - Contains CUDA compression code								//
//																				//
//  This is a GPU implementation of a compressor for the PVR format as			//
//  described by Simon Fenney in his paper Texture Compression using			//
//  Low-Frequency Signal Modulation.											//
//																				//
//////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////////
//									File Includes	 							//
//////////////////////////////////////////////////////////////////////////////////
#include <cuda.h>
#include <cuda_runtime_api.h>
#include <cutil_inline.h>
#include <stdio.h>
#include <limits.h>
//////////////////////////////////////////////////////////////////////////////////
//									File Includes	 							//
//////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////////
//							   Preprocessor Definitions	 						//
//////////////////////////////////////////////////////////////////////////////////
#define USE_PIXEL_UPDATE		// Method to calculate 
#define USE_SVD					// Use SVD or MP to perform optimization
#define GET_RMS_ERROR			// Whether or not to return RMS error
//#define DECOMPRESS_PVR		// Whether the final image is decompressed
#define TWO_BY_TWO				// Size of the optimization block (2x2 or 1x1)
#define BLOCK_WIDTH				16
#define BLOCK_HEIGHT			16
#define FILTER_LENGTH			2
#define HALF_FILTER_LENGTH		FILTER_LENGTH * 0.5f
#define RED_SHIFT				16
#define GREEN_SHIFT				8
#define	BLUE_SHIFT				0
#define RSHIFT_16BPP			3
#define LSHIFT_16BPP			5
#define ALPHA_MASK				0xFF000000
#define RED_MASK				0x00FF0000
#define GREEN_MASK				0x0000FF00
#define BLUE_MASK				0x000000FF
#define ONE_EIGHTH				0.125f
#define ONE_FOURTH				0.25f
#define THREE_EIGHTHS			0.375f
#define FIVE_EIGHTHS			0.625f
#define NUM_OPTIMIZATION_PASSES	20
#define EPSILON					10e-10
//////////////////////////////////////////////////////////////////////////////////
//							   Preprocessor Definitions	 						//
//////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////////
//								  Macro Definitions		 						//
//////////////////////////////////////////////////////////////////////////////////
#define SQR(x)					((x) * (x))
#define CLAMP(x, a, b)			(min((max((x), (a))), (b)))
#define MAKE_RED_PIXEL(x)		(((x) & RED_MASK)>>RED_SHIFT)
#define MAKE_GREEN_PIXEL(x)		(((x) & GREEN_MASK)>>GREEN_SHIFT)
#define MAKE_BLUE_PIXEL(x)		(((x) & BLUE_MASK)>>BLUE_SHIFT)
#define MAKE_ARGB(r, g, b)		(ALPHA_MASK | ((r)<<RED_SHIFT) | ((g)<<GREEN_SHIFT) | ((b)<<BLUE_SHIFT))
//////////////////////////////////////////////////////////////////////////////////
//								  Macro Definitions		 						//
//////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////////
//									   Misc.									//
//////////////////////////////////////////////////////////////////////////////////
// CUDA texture declarations
texture<unsigned int, cudaTextureType2D, cudaReadModeElementType> texRef;
texture<unsigned int, cudaTextureType2D, cudaReadModeElementType> origRef;
texture<int, cudaTextureType2D, cudaReadModeElementType> redCurrentRef;
texture<int, cudaTextureType2D, cudaReadModeElementType> greenCurrentRef;
texture<int, cudaTextureType2D, cudaReadModeElementType> blueCurrentRef;
texture<float, cudaTextureType2D, cudaReadModeElementType> modRef;
texture<float, cudaTextureType1D, cudaReadModeElementType> filterRef;
#ifdef TWO_BY_TWO
texture<float, cudaTextureType1D, cudaReadModeElementType> svdMatRefTL;
texture<float, cudaTextureType1D, cudaReadModeElementType> svdMatRefTR;
texture<float, cudaTextureType1D, cudaReadModeElementType> svdMatRefBL;
texture<float, cudaTextureType1D, cudaReadModeElementType> svdMatRefBR;
#else
texture<float, cudaTextureType1D, cudaReadModeElementType> svdMatRef;
#endif // TWO_BY_TWO

// various wavelet filters. remember to update the FILTER_LENGTH macro if you
// choose to use a different one

// bior wavelet filter
//static const float wavelet_filter[] = {
//	-0.00488281250000000000,
//	0.00976562500000000000,
//	0.03320312500000000000,
//	-0.07617187500000000000,
//	-0.12011718750000000000,
//	0.31640625000000000000,
//	0.68359375000000000000,
//	0.31640625000000000000,
//	-0.12011718750000000000,
//	-0.07617187500000000000,
//	0.03320312500000000000,
//	0.00976562500000000000,
//	-0.00488281250000000000,
//};

// 2nd bior wavelet filter
//static const float wavelet_filter[] = {
//-0.0625,
//0.06250,
//0.50000,
//0.50000,
//0.06250,
//-0.0625
//};

// my wavelet filter
//static const float wavelet_filter[] = {
//-0.2500,
//0.50000,
//0.50000,
//0.50000,
//-0.2500
//};

// daubechies wavelet filter
//static const float wavelet_filter[] = {
//	0.02692517479416041400,
//	0.17241715192471294000,
//	0.42767453217028290000,
//	0.46477285717277800000,
//	0.09418477475112015100,
//	-0.20737588089628295000,
//	-0.06847677451090331000,
//	0.10503417113713563000,
//	0.02172633772990401800,
//	-0.04782363205881859400,
//	0.00017744640673182261,
//	0.01581208292613723100,
//	-0.00333981011324138060,
//	-0.00302748028715121120,
//	0.00130648364017893680,
//	0.00016290733600968354,
//	-0.00017816487954739422,
//	0.00002782275679290904
//};

// simple wavelet filter
static float wavelet_filter[] = {0.5f, 0.5f};

// weight matrix to be combined with modulation data in the optimization step
#define ROW(n)	1*n,	2*n,	3*n,	4*n,	3*n,	2*n,	1*n
#ifdef TWO_BY_TWO
#define USE_SVD
#define SVD_FACTOR				ONE_EIGHTH
#define SVD_OFFSET				4
#define SVD_MAT_WIDTH			8
#define SVD_MAT_HEIGHT			121
#define SVD_DIM					11
static const float MwTL[SVD_MAT_HEIGHT] = 
{
	ROW(0.0625f), 0, 0, 0, 0,	// 1/16
	ROW(0.125f), 0, 0, 0, 0,	// 2/16
	ROW(0.1875f), 0, 0, 0, 0,	// 3/16
	ROW(0.25f), 0, 0, 0, 0,		// 4/16
	ROW(0.1875f), 0, 0, 0, 0,	// 3/16
	ROW(0.125f), 0, 0, 0, 0,	// 2/16
	ROW(0.0625f), 0, 0, 0, 0,	// 1/16
	ROW(0), 0, 0, 0, 0,
	ROW(0), 0, 0, 0, 0,
	ROW(0), 0, 0, 0, 0,
	ROW(0), 0, 0, 0, 0
};
static const float MwTR[SVD_MAT_HEIGHT] = 
{
	0, 0, 0, 0, ROW(0.0625f),	// 1/16
	0, 0, 0, 0, ROW(0.125f),	// 2/16
	0, 0, 0, 0, ROW(0.1875f),	// 3/16
	0, 0, 0, 0, ROW(0.25f),		// 4/16
	0, 0, 0, 0, ROW(0.1875f),	// 3/16
	0, 0, 0, 0, ROW(0.125f),	// 2/16
	0, 0, 0, 0, ROW(0.0625f),	// 1/16
	ROW(0), 0, 0, 0, 0,
	ROW(0), 0, 0, 0, 0,
	ROW(0), 0, 0, 0, 0,
	ROW(0), 0, 0, 0, 0
};
static const float MwBL[SVD_MAT_HEIGHT] = 
{
	ROW(0), 0, 0, 0, 0,
	ROW(0), 0, 0, 0, 0,
	ROW(0), 0, 0, 0, 0,
	ROW(0), 0, 0, 0, 0,
	ROW(0.0625f), 0, 0, 0, 0,	// 1/16
	ROW(0.125f), 0, 0, 0, 0,	// 2/16
	ROW(0.1875f), 0, 0, 0, 0,	// 3/16
	ROW(0.25f), 0, 0, 0, 0,		// 4/16
	ROW(0.1875f), 0, 0, 0, 0,	// 3/16
	ROW(0.125f), 0, 0, 0, 0,	// 2/16
	ROW(0.0625f), 0, 0, 0, 0,	// 1/16
};
static const float MwBR[SVD_MAT_HEIGHT] = 
{
	ROW(0), 0, 0, 0, 0,
	ROW(0), 0, 0, 0, 0,
	ROW(0), 0, 0, 0, 0,
	ROW(0), 0, 0, 0, 0,
	0, 0, 0, 0, ROW(0.0625f),	// 1/16
	0, 0, 0, 0, ROW(0.125f),	// 2/16
	0, 0, 0, 0, ROW(0.1875f),	// 3/16
	0, 0, 0, 0, ROW(0.25f),		// 4/16
	0, 0, 0, 0, ROW(0.1875f),	// 3/16
	0, 0, 0, 0, ROW(0.125f),	// 2/16
	0, 0, 0, 0, ROW(0.0625f)	// 1/16
};
#else
#define SVD_FACTOR				ONE_FOURTH
#define SVD_OFFSET				4
#define SVD_MAT_WIDTH			2
#define SVD_MAT_HEIGHT			49
#define SVD_DIM					7
static const float Mw[SVD_MAT_HEIGHT] = 
{
	ROW(0.0625f),	// 1/16
	ROW(0.125f),	// 2/16
	ROW(0.1875f),	// 3/16
	ROW(0.25f),		// 4/16
	ROW(0.1875f),	// 3/16
	ROW(0.125f),	// 2/16
	ROW(0.0625f)	// 1/16
};
#endif // TWO_BY_TWO
// pointer swap method
static inline void swap(void **x, void **y) {
	void *t = *x;
	*x = *y;
	*y = t;
}

//////////////////////////////////////////////////////////////////////////////////
//	Calculate the root mean squared error of the image that is decoded from the	//
//	compressed data. This kernel calculates the total error for the red, green,	//
//	and blue channels of a single pixel and stores it in an array.				//
//																				//
//	orig		Pointer to the pixel data for the original image				//
//	a			Pointer to the upscaled pixel data for the final A image		//
//	b			Pointer to the upscaled pixel data for the final B image		//
//	out			Pointer to the array where the total error for the current 		//
//				pixel is stored													//
//	width		The width, in pixels, of the original image						//
//																				//
//////////////////////////////////////////////////////////////////////////////////
#ifdef GET_RMS_ERROR
__global__ void rms_error(unsigned int *orig, unsigned int *a, unsigned int *b,
						  unsigned int *out,  int width) {
	unsigned int y = blockDim.y * blockIdx.y + threadIdx.y;
	unsigned int x = blockDim.x * blockIdx.x + threadIdx.x;
	int idx = (width * y) + x;

	unsigned int opixel = orig[idx];
	unsigned int apixel = a[idx];
	unsigned int bpixel = b[idx];
	float modbit = tex2D(modRef, x, y);
	float r_modbit = 1.0f - modbit;
	int c_red = ((float)MAKE_RED_PIXEL(bpixel)*modbit + 
		(float)MAKE_RED_PIXEL(apixel)*r_modbit);
	int c_green = ((float)MAKE_GREEN_PIXEL(bpixel)*modbit + 
		(float)MAKE_GREEN_PIXEL(apixel)*r_modbit);
	int c_blue = ((float)MAKE_BLUE_PIXEL(bpixel)*modbit + 
		(float)MAKE_BLUE_PIXEL(apixel)*r_modbit);

	out[idx] = SQR(((float)MAKE_RED_PIXEL(opixel) - c_red)) + 
			   SQR(((float)MAKE_GREEN_PIXEL(opixel) - c_green)) + 
			   SQR(((float)MAKE_BLUE_PIXEL(opixel) - c_blue));
}
#endif // GET_RMS_ERROR

#ifdef DECOMPRESS_PVR
__global__ void decompress(unsigned int *a, unsigned int *b, 
						   unsigned int *out, float *mod, int width) {
	unsigned int y = blockDim.y * blockIdx.y + threadIdx.y;
	unsigned int x = blockDim.x * blockIdx.x + threadIdx.x;
	int idx = (width * y) + x;

	unsigned int apixel = a[idx];
	int a_red = MAKE_RED_PIXEL(apixel);
	int a_green = MAKE_GREEN_PIXEL(apixel);
	int a_blue = MAKE_BLUE_PIXEL(apixel);

	unsigned int bpixel = b[idx];
	int b_red = MAKE_RED_PIXEL(bpixel);
	int b_green = MAKE_GREEN_PIXEL(bpixel);
	int b_blue = MAKE_BLUE_PIXEL(bpixel);

	float modbit = mod[idx];
	float r_modbit = 1.0f - modbit;

	int c_red = ((float)b_red*modbit + (float)a_red*r_modbit);
	int c_green = ((float)b_green*modbit + (float)a_green*r_modbit);
	int c_blue = ((float)b_blue*modbit + (float)a_blue*r_modbit);

	out[idx] = MAKE_ARGB(c_red, c_green, c_blue);

}
#endif // DECOMPRESS_PVR

__global__ void encode_texture(const unsigned int* A, unsigned int *B, float *mod, 
							   int width, bool hasAlpha, unsigned int *out, 
							   int outWidth, int outHeight) {
	unsigned int y = blockDim.y * blockIdx.y + threadIdx.y;
	unsigned int x = blockDim.x * blockIdx.x + threadIdx.x;
	unsigned outy = y, outx = x;
	//int idx = (width * CLAMP(y, 15, outHeight-15)) + (CLAMP(x, 15, outWidth-15));
	/*if (x < 10) {
		outx = 15;
		outy = 15;
	} else if (x > outWidth-10) {
		outx = 80;
		outy = 45;
	} 
	if (y < 10) {
		outy = 15;
		outx = 15;
	} else if (y > outHeight-10) {
		outy = 45;
		outx = 80;
	}*/ 
	int idx = (width*outy) + outx;

	// Interleave lower 16 bits of x and y, so the bits of y
	// are in the even positions and bits from x in the odd;
	// outIdx is the resulting 32-bit Morton Number. 
	//unsigned int a = (2*x >= outWidth) ? 2*x - outWidth : 2*x;
	unsigned int a = 2*x;
	unsigned int b = y;

	a = (a | (a << 8)) & 0x00FF00FF;
	a = (a | (a << 4)) & 0x0F0F0F0F;
	a = (a | (a << 2)) & 0x33333333;
	a = (a | (a << 1)) & 0x55555555;

	b = (b | (b << 8)) & 0x00FF00FF;
	b = (b | (b << 4)) & 0x0F0F0F0F;
	b = (b | (b << 2)) & 0x33333333;
	b = (b | (b << 1)) & 0x55555555;

	// multiply by two because we encode 64-bit blocks not 32-bit blocks
	//int outIdx = 2 * CLAMP((a | (b << 1)), 0, outWidth*outHeight);
	// TODO: Fix this, currently only works if dimnesion is 2x as big
	int outIdx = (a | (b << 1));
	//outIdx += (2*x >= outWidth) ? (outWidth*outHeight) : 0;

	unsigned int apixel = A[idx];
	unsigned int bpixel = B[idx];
	unsigned int color = 0, modbits = 0;

	if (hasAlpha == false) {
		// encode color B
		color = 1;
		color = ((color<<LSHIFT_16BPP) | (MAKE_BLUE_PIXEL(bpixel)>>RSHIFT_16BPP));
		color = ((color<<LSHIFT_16BPP) | (MAKE_GREEN_PIXEL(bpixel)>>RSHIFT_16BPP));
		color = ((color<<LSHIFT_16BPP) | (MAKE_RED_PIXEL(bpixel)>>RSHIFT_16BPP));

		// encode color A
		color = ((color<<1) | 1);
		color = ((color<<LSHIFT_16BPP) | (MAKE_BLUE_PIXEL(apixel)>>RSHIFT_16BPP));
		color = ((color<<LSHIFT_16BPP) | (MAKE_GREEN_PIXEL(apixel)>>RSHIFT_16BPP));
		color = ((color<<(LSHIFT_16BPP-1)) | (MAKE_RED_PIXEL(apixel)>>(RSHIFT_16BPP+1)));

		// modulation mode 0
		color = ((color<<1) | 0);

		// encode modulation bits
		int i, j;
		float mbit;
		for (i = 3; i >= 0; i--) {
			for (j = 3; j >= 0; j--) {
				mbit = tex2D(modRef, j + (4*x), i + (4*y));
				if (mbit < 0.05f) {
					modbits = (modbits<<2);
				} else if (mbit < 0.5f) {
					modbits = ((modbits<<2) | (0x1));
				} else if (mbit < 0.95f) {
					modbits = ((modbits<<2) | (0x2));
				} else {
					modbits = ((modbits<<2) | (0x3));
				}
			}
		}

		out[outIdx] = modbits;
		out[outIdx+1] = color;
	} else {
		// TODO: implement this
	}
}
//////////////////////////////////////////////////////////////////////////////////
//									   Misc.									//
//////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////////
//								  Optimization									//
//////////////////////////////////////////////////////////////////////////////////
__device__ float pythag(float a, float b) {
	float absa = fabsf(a);
	float absb = fabsf(b);
	if (absa > absb) return absa*sqrtf(1.0f+SQR((absb/absa)));
	else return (absb == 0.0f ? 0.0f : absb*sqrtf(1.0f+SQR(absa/absb)));
}

__device__ float sign(float a, float b) {
	return b >= 0.0f ? ( a >= 0.0f ? a : -a) : ( a >= 0.0f ? -a : a);
}

#ifndef TWO_BY_TWO
//////////////////////////////////////////////////////////////////////////////////
//	Calculate the Moore-Penrose pseudo-inverse of the weight matrix and use		//
//	it to compute the updated pair of A and B representative pixels.			//
//																				//
//	candidateA	Pointer to the output location of the A image					//
//	candidateB	Pointer to the output location of the B image					//
//	oldA		Pointer to the A image computed in the previous iteration		//
//	oldB		Pointer to the B image computed in the previous iteration		//
//	width		The width, in pixels, of the image to be compressed				//
//	err			Pointer to memory where errors are reported						//
//																				//
//////////////////////////////////////////////////////////////////////////////////
__global__ void moore_penrose_optimize(unsigned int *candidateA, unsigned int *candidateB,
									   unsigned int *oldA, unsigned int *oldB, 
									   int width, int height, int *err) {
	int thready = blockDim.y * blockIdx.y + threadIdx.y;
	int threadx = blockDim.x * blockIdx.x + threadIdx.x;
	const int  tIdx = (width*thready) + threadx;
	int red,red0=0,red1=0,green,green0=0,green1=0,blue,blue0=0,blue1=0;
	int x_offset = SVD_OFFSET*threadx-1, y_offset = SVD_OFFSET*thready-1;
	int x, y, pixelx, pixely, index;
	float A00=0,A01=0,A10=0,A11=0,InverseA00,InverseA01,InverseA10,InverseA11;
	float a0, a1, modbit, dist, det;

	// Fetch the weight matrix for the 7x7 optimization window of the current
	// pair of representative values
	for (y = 0; y < SVD_DIM; y++) {
		for (x = 0; x < SVD_DIM; x++) {
			index = y*SVD_DIM + x;
			pixelx = x_offset + x;
			pixely = y_offset + y;

			// fetch all the necessary values
			modbit = tex2D(modRef, CLAMP(pixelx, 0, width - 2), CLAMP(pixely, 0, height-2));
			red = tex2D(redCurrentRef, CLAMP(pixelx, 0, width-2), CLAMP(pixely, 0, height-2));
			green = tex2D(greenCurrentRef, CLAMP(pixelx, 0, width-2), CLAMP(pixely, 0, height-2));
			blue = tex2D(blueCurrentRef, CLAMP(pixelx, 0, width-2), CLAMP(pixely, 0, height-2));
			dist = tex1Dfetch(svdMatRef, index);
			a0 = dist*(1.0f - modbit);
			a1 = dist*modbit;	

			// A = TransposeA * A
			A00 += a0*a0;
			A01 += a0*a1;
			A10 += a1*a0;
			A11 += a1*a1;

			// colors = TransposeA * color
			red0 += a0*red;
			red1 += a1*red;

			green0 += a0*green;
			green1 += a1*green;

			blue0 += a0*blue;
			blue1 += a1*blue;
		}
	}

	// since we only have a 2x2 matrix, manually calculate inverse
	det = 1.0f / (A00*A11 - A01*A10);
	InverseA00 = det * A11;
	InverseA01 = -1 * det * A01;
	InverseA10 = -1 * det * A10;
	InverseA11 = det * A00;

	// Calculate new representative colors by using a "fix the error" approach.
	// We multiply the pseudo-inverse weight matrix by the difference between the
	// current and original images to get an "update" that is applied to each
	// representative value. We also limit how big the update can be to avoid
	// colors flying out of bounds.
	unsigned int oldColor;

	// A representative
	oldColor = oldA[tIdx];
	red = InverseA00*(float)red0 + InverseA01*(float)red1;
	green = InverseA00*(float)green0 + InverseA01*(float)green1;
	blue = InverseA00*(float)blue0 + InverseA01*(float)blue1;
#ifdef USE_PIXEL_UPDATE
	candidateA[tIdx] = MAKE_ARGB(
		CLAMP((int)MAKE_RED_PIXEL(oldColor) + CLAMP((int)red, -16, 16), 0, 255),
		CLAMP((int)MAKE_GREEN_PIXEL(oldColor) + CLAMP((int)green, -16, 16), 0, 255),
		CLAMP((int)MAKE_BLUE_PIXEL(oldColor) + CLAMP((int)blue, -16, 16), 0, 255));
#else
	candidateA[tIdx] = MAKE_ARGB(CLAMP(red, 0, 255), CLAMP(green, 0, 255), CLAMP(blue, 0, 255));
#endif // USE_PIXEL_UPDATE
	oldA[tIdx] = candidateA[tIdx];

	// B representative
	oldColor = oldB[tIdx];
	red = InverseA10*(float)red0 + InverseA11*(float)red1;
	green = InverseA10*(float)green0 + InverseA11*(float)green1;
	blue = InverseA10*(float)blue0 + InverseA11*(float)blue1;
#ifdef USE_PIXEL_UPDATE
	candidateB[tIdx] = MAKE_ARGB(
		CLAMP((int)MAKE_RED_PIXEL(oldColor) + CLAMP((int)red, -16, 16), 0, 255),
		CLAMP((int)MAKE_GREEN_PIXEL(oldColor) + CLAMP((int)green, -16, 16), 0, 255),
		CLAMP((int)MAKE_BLUE_PIXEL(oldColor) + CLAMP((int)blue, -16, 16), 0, 255));
#else
	candidateB[tIdx] = MAKE_ARGB(CLAMP(red, 0, 255), CLAMP(green, 0, 255), CLAMP(blue, 0, 255));
#endif // USE_PIXEL_UPDATE
	oldB[tIdx] = candidateB[tIdx];
}
#endif // TWO_BY_TWO
//////////////////////////////////////////////////////////////////////////////////
//	Calculates the singular value decomposition of the weight matrix and uses	//
//	it to compute teh pseudo-inverse. Then multiplies the pseudo-inverse by		//
//	the current color values to get the new values.								//
//	DOCUMENTATION: Numerical Recipes in C, 2nd Edition, chapter 2.6, 11.2,		//
//	& 11.3																		//
//																				//
//	candidateA	Pointer to the output location of the A image					//
//	candidateB	Pointer to the output location of the B image					//
//	oldA		Pointer to the A image computed in the previous iteration		//
//	oldB		Pointer to the B image computed in the previous iteration		//
//	width		The width, in pixels, of the image to be compressed				//
//	err			Pointer to memory where errors are reported						//
//																				//
//////////////////////////////////////////////////////////////////////////////////
__global__ void svd_optimize(unsigned int *candidateA, unsigned int *candidateB,
									unsigned int *oldA, unsigned int *oldB, 
									int width,int height,  int *err) {
	int thready = 2 * (blockDim.y * blockIdx.y + threadIdx.y);
	int threadx = 2 * (blockDim.x * blockIdx.x + threadIdx.x);
	const int maxits=30;
	int tIdx = (width*thready) + threadx;
	int i,its,j,jj,k,l=0,nm=0,red[SVD_MAT_HEIGHT],green[SVD_MAT_HEIGHT],
		blue[SVD_MAT_HEIGHT];
	bool flag;
	float anorm,c,f,g,h,s,scale,w,x,y,z,rv[SVD_MAT_WIDTH]; 
	// column major order for SVD
	float A[SVD_MAT_HEIGHT][SVD_MAT_WIDTH], Eigenvalues[SVD_MAT_WIDTH], 
		Eigenbasis[SVD_MAT_WIDTH][SVD_MAT_WIDTH], 
		InverseA[SVD_MAT_WIDTH][SVD_MAT_HEIGHT];

#ifdef USE_SHARED_MEM
	// WARNING - DOES NOT COMPILE. TOO MUCH SHARED MEMORY
	__shared__ int sRed[4*(BLOCK_HEIGHT+1)][4*(BLOCK_WIDTH+1)], 
				   sGreen[4*(BLOCK_HEIGHT+1)][4*(BLOCK_WIDTH+1)], 
				   sBlue[4*(BLOCK_HEIGHT+1)][4*(BLOCK_WIDTH+1)];
	__shared__ float sMod[4*(BLOCK_HEIGHT+1)][4*(BLOCK_WIDTH+1)], 
					 sDist[SVD_MAT_HEIGHT];

	bool maxX = false, maxY = false;
	int shrI, shrJ, shrX, shrY, shrtx = 4*threadIdx.x, shrty = 4*threadIdx.y;
	for (shrY = 0, shrI = (SVD_OFFSET*thready) - 1; shrY < 4; shrY++, shrI++) {
		for (shrX = 0, shrJ = (SVD_OFFSET*threadx) - 1; shrX < 4; shrX++, shrJ++) {
			sRed[shrty + shrY][shrtx + shrX] = tex2D(redCurrentRef, shrJ, shrI);
			sGreen[shrty + shrY][shrtx + shrX] = tex2D(greenCurrentRef, shrJ, shrI);
			sBlue[shrty + shrY][shrtx + shrX] = tex2D(blueCurrentRef, shrJ, shrI);
			sMod[shrty + shrY][shrtx + shrX] = tex2D(modRef, shrJ, shrI);
		}
	}

	// check if we are at the x border
	if (threadIdx.x == blockDim.x - 1) {
		maxX = true;
		shrtx = 4*(threadIdx.x+1);
		shrty = 4*threadIdx.y;
		for (shrY = 0, shrI = (SVD_OFFSET*thready) - 1; shrY < 4; shrY++, shrI++) {
			for (shrX = 0, shrJ = (SVD_OFFSET*(threadx+1)) - 1; shrX < 4; shrX++, shrJ++) {
				sRed[shrty + shrY][shrtx + shrX] = tex2D(redCurrentRef, shrJ, shrI);
				sGreen[shrty + shrY][shrtx + shrX] = tex2D(greenCurrentRef, shrJ, shrI);
				sBlue[shrty + shrY][shrtx + shrX] = tex2D(blueCurrentRef, shrJ, shrI);
				sMod[shrty + shrY][shrtx + shrX] = tex2D(modRef, shrJ, shrI);
			}
		}
	}

	// check if we are at the y border
	if (threadIdx.y == blockDim.y - 1) {
		maxY = true;
		shrtx = 4*threadIdx.x;
		shrty = 4*(threadIdx.y+1);
		for (shrY = 0, shrI = (SVD_OFFSET*(thready+1)) - 1; shrY < 4; shrY++, shrI++) {
			for (shrX = 0, shrJ = (SVD_OFFSET*threadx) - 1; shrX < 4; shrX++, shrJ++) {
				sRed[shrty + shrY][shrtx + shrX] = tex2D(redCurrentRef, shrJ, shrI);
				sGreen[shrty + shrY][shrtx + shrX] = tex2D(greenCurrentRef, shrJ, shrI);
				sBlue[shrty + shrY][shrtx + shrX] = tex2D(blueCurrentRef, shrJ, shrI);
				sMod[shrty + shrY][shrtx + shrX] = tex2D(modRef, shrJ, shrI);
			}
		}
	}

	// check if we are at the corner
	if (maxX && maxY) {
		shrtx = 4*(threadIdx.x+1);
		shrty = 4*(threadIdx.y+1);
		for (shrY = 0, shrI = (SVD_OFFSET*(thready+1)) - 1; shrY < 4; shrY++, shrI++) {
			for (shrX = 0, shrJ = (SVD_OFFSET*(threadx+1)) - 1; shrX < 4; shrX++, shrJ++) {
				sRed[shrty + shrY][shrtx + shrX] = tex2D(redCurrentRef, shrJ, shrI);
				sGreen[shrty + shrY][shrtx + shrX] = tex2D(greenCurrentRef, shrJ, shrI);
				sBlue[shrty + shrY][shrtx + shrX] = tex2D(blueCurrentRef, shrJ, shrI);
				sMod[shrty + shrY][shrtx + shrX] = tex2D(modRef, shrJ, shrI);
			}
		}
	}

	// check if we have to fetch the distance matrix
	if (tIdx < 49) {
		sDist[tIdx] = tex1Dfetch(svdMatRef, tIdx);
	}
	__syncthreads();
	
	// now get matrix A
	int pixelx, pixely;
	float dist, modbit;
	shrtx = 4*threadIdx.x;
	shrty = 4*threadIdx.y;
	for (shrJ = 0; shrJ < SVD_MAT_HEIGHT; shrJ++) {
		// j is the index of the pixel we want
		pixelx = (shrJ % SVD_DIM) + shrtx;
		pixely = (shrJ / SVD_DIM) + shrty;
		modbit = sMod[pixely][pixelx];
		red[j] = sRed[pixely][pixelx];
		green[j] = sGreen[pixely][pixelx];
		blue[j] = sBlue[pixely][pixelx];
		dist = sDist[j];
		A[j][0] = dist*(1.0f - modbit);
		A[j][1] = dist*modbit;	
	}
#else	// Not using shared memory
	// Fetch the weight matrix for the optimization window of the current
	// pair of representative pixels
	int index, pixelx, pixely;
	float dist, modbit;
	for (j = 0; j < SVD_DIM; j++) {
		for (i = 0; i < SVD_DIM; i++) {
			// j is the index of the pixel we want
			index = j*SVD_DIM + i;
			pixelx = SVD_OFFSET*threadx-1 + i;
			pixely = SVD_OFFSET*thready-1 + j;

			// fetch all the necessary values
			modbit = tex2D(modRef, CLAMP(pixelx, 0, width), CLAMP(pixely, 0, height));
			red[index] = tex2D(redCurrentRef, CLAMP(pixelx, 0, width), CLAMP(pixely, 0, height));
			green[index] = tex2D(greenCurrentRef, CLAMP(pixelx, 0, width), CLAMP(pixely, 0, height));
			blue[index] = tex2D(blueCurrentRef, CLAMP(pixelx, 0, width), CLAMP(pixely, 0, height));
#ifdef TWO_BY_TWO
			dist = tex1Dfetch(svdMatRefTL, index);
			A[index][0] = dist*(1.0f - modbit);
			A[index][1] = dist*modbit;
			dist = tex1Dfetch(svdMatRefTR, index);
			A[index][2] = dist*(1.0f - modbit);
			A[index][3] = dist*modbit;
			dist = tex1Dfetch(svdMatRefBL, index);
			A[index][4] = dist*(1.0f - modbit);
			A[index][5] = dist*modbit;
			dist = tex1Dfetch(svdMatRefBR, index);
			A[index][6] = dist*(1.0f - modbit);
			A[index][7] = dist*modbit;
#else
			dist = tex1Dfetch(svdMatRef, index);
			A[index][0] = dist*(1.0f - modbit);
			A[index][1] = dist*modbit;	
#endif //TWO_BY_TWO
		}
	}
#endif	// USE_SHARED_MEM

	// Householder reduction to bidiagonal form
	g = scale = anorm = 0.0f;
	for (i=0; i<SVD_MAT_WIDTH; i++) {
		l = i+1;
		rv[i] = scale*g;
		g = s = scale = 0.0f;
		if (i<SVD_MAT_HEIGHT) {
			for (k=i; k<SVD_MAT_HEIGHT; k++) scale+= fabsf(A[k][i]);
			if (scale!=0.0f) {
				for (k=i; k<SVD_MAT_HEIGHT; k++) {
					A[k][i]/= scale;
					s+= SQR(A[k][i]);
				}
				f = A[i][i];
				g = -sign(sqrtf(s),f);
				h = f*g-s;
				A[i][i] = f-g;
				for (j=l; j<SVD_MAT_WIDTH; j++) {
				  for (s=0.0f,k=i; k<SVD_MAT_HEIGHT; k++) s+= A[k][i]*A[k][j];
					f = s/h;
					for (k=i; k<SVD_MAT_HEIGHT; k++) A[k][j]+= f*A[k][i];
				}
				for (k=i; k<SVD_MAT_HEIGHT; k++) A[k][i]*= scale;
			}
		}
		Eigenvalues[i] = scale*g;
		g = s = scale = 0.0f;
		if (i<SVD_MAT_HEIGHT && i!=SVD_MAT_WIDTH-1) {
			for (k=l; k<SVD_MAT_WIDTH; k++) scale+= fabsf(A[i][k]);
			if (scale!=0.0f)  {
				for(k=l; k<SVD_MAT_WIDTH; k++) {
					A[i][k]/= scale;
					s+= SQR(A[i][k]);
				}
				f = A[i][l];
				g = -sign(sqrtf(s),f);
				h = f*g-s;
				A[i][l] = f-g;
				for (k=l; k<SVD_MAT_WIDTH; k++) rv[k] = A[i][k]/h;
				for (j=l; j<SVD_MAT_HEIGHT; j++) {
				  for(s=0.0f,k=l; k<SVD_MAT_WIDTH; k++) s+= A[j][k]*A[i][k];
					for(k=l; k<SVD_MAT_WIDTH; k++) A[j][k]+= s*rv[k];
				}
				for(k=l; k<SVD_MAT_WIDTH; k++) A[i][k]*= scale;
			}
		}
		anorm = fmaxf(anorm,(fabsf(Eigenvalues[i])+fabsf(rv[i])));
	}
	
	// Accumulate right-hand side updates
	for(i=SVD_MAT_WIDTH-1; i>=0; i--) {
		if (i<SVD_MAT_WIDTH-1) {
		  if (g!=0.0f) {
				for (j=l; j<SVD_MAT_WIDTH; j++) Eigenbasis[j][i] = (A[i][j]/A[i][l])/g;
				for (j=l; j<SVD_MAT_WIDTH; j++) {
				  for(s=0.0f,k=l; k<SVD_MAT_WIDTH; k++) s+= A[i][k]*Eigenbasis[k][j];
					for(k=l; k<SVD_MAT_WIDTH; k++) Eigenbasis[k][j]+= s*Eigenbasis[k][i];
				}
			}
			for (j=l; j<SVD_MAT_WIDTH; j++) Eigenbasis[i][j] = Eigenbasis[j][i] = 0.0f;
		}
		Eigenbasis[i][i] = 1.0f;
		g = rv[i];
		l = i;
	}
		
	// Accumulate left-hand side updates
	for (i=min(SVD_MAT_HEIGHT-1,SVD_MAT_WIDTH-1); i>=0; i--) {
		l = i+1;
		g = Eigenvalues[i];
		for(j=l;j<SVD_MAT_WIDTH;j++) A[i][j] = 0.0f;
		if (g!=0.0f) {
		  g = 1.0f/g;
			for (j=l; j<SVD_MAT_WIDTH; j++) {
			  for (s=0.0f,k=l; k<SVD_MAT_HEIGHT; k++) s+= A[k][i]*A[k][j];
				f = (s/A[i][i])*g;
				for (k=i; k<SVD_MAT_HEIGHT; k++) A[k][j]+= f*A[k][i];
			}
			for (j=i; j<SVD_MAT_HEIGHT; j++) A[j][i]*= g;
		}
		else for (j=i; j<SVD_MAT_HEIGHT; j++) A[j][i] = 0.0f;
		A[i][i]+= 1.0f;
	}
	
	// diagonalization of the bidiagonal form: loop over singular values, and
	// over allowed iterations
	for (k=SVD_MAT_WIDTH-1; k>=0; k--) {
		for (its=0; its<maxits; its++) {
			flag = true;
			for (l=k; l>=0; l--) {			// test for splitting
				nm = l-1;
				if ((float)(fabsf(rv[l])+anorm) == anorm) {
					flag =  false;
					break;
				}
				if ((float)(fabsf(Eigenvalues[nm])+anorm) == anorm) break;
			}
			if (flag) {
			  c = 0.0f;
			  s = 1.0f;
				for (i=l; i<=k; i++) {
					f = s*rv[i];
					rv[i]*= c;
					if ((float)(fabsf(f)+anorm) == anorm) break;
					g = Eigenvalues[i];
					h = pythag(f,g);
					Eigenvalues[i] = h;
					h = 1.0f/h;
					c = g*h;
					s = -f*h;
					for (j=0; j<SVD_MAT_HEIGHT; j++) {
						y = A[j][nm];
						z = A[j][i];
						A[j][nm] = y*c+z*s;
						A[j][i]  = z*c-y*s;
					}
				}
			}
			z = Eigenvalues[k];
			if (l==k) {					// convergence
				if(z<0.0f) {			// singular value is made non-negative
					Eigenvalues[k] = -z;
					for(j=0; j<SVD_MAT_WIDTH; j++) Eigenbasis[j][k] *= -1.0f;
				}
				break;
			}
			
			if(its>=maxits) { *err = 1; } // error check
			x = Eigenvalues[l];			// shift from bottom 2-by-2 minor
			nm = k-1;
			y = Eigenvalues[nm];
			g = rv[nm];
			h = rv[k];
			f = ((y-z)*(y+z)+(g-h)*(g+h))/(2.0f*h*y);
			g = pythag(f,1.0);
			f = ((x-z)*(x+z)+h*((y/(f+sign(g,f)))-h))/x;
			c = s = 1.0f;				// next QR transformation
			for (j=l; j<=nm; j++) {
				i = j+1;
				g = rv[i];
				y = Eigenvalues[i];
				h = s*g;
				g = c*g;
				z = pythag(f,h);
				rv[j] = z;
				c = f/z;
				s = h/z;
				f = x*c+g*s;
				g = g*c-x*s;
				h = y*s;
				y*= c;
				for(jj=0;jj<SVD_MAT_WIDTH;jj++) {
					x = Eigenbasis[jj][j];
					z = Eigenbasis[jj][i];
					Eigenbasis[jj][j] = x*c+z*s;
					Eigenbasis[jj][i] = z*c-x*s;
				}
				z = pythag(f,h);
				Eigenvalues[j] = z;		// rotation can be arbitrary if z = 0
				if (z!=0.0f) {
				  z = 1.0f / z;
					c = f*z;
					s = h*z;
				}
				f = c*g+s*y;
				x = c*y-s*g;
				for (jj=0; jj<SVD_MAT_HEIGHT; jj++) {
					y = A[jj][j];
					z = A[jj][i];
					A[jj][j] = y*c+z*s;
					A[jj][i] = z*c-y*s;
				}
			}
			rv[l] = 0.0f;
			rv[k] = f;
			Eigenvalues[k] = x;
		}
	}

	// begin constructing inverse matrix
	// inverse(A) = V * inverse(W) * transpose(A)

	// Step 1: transpose(A) = inverse(W) * transpose(A)
	for (i = 0; i < SVD_MAT_WIDTH; i++) {
		w = 1.0f / Eigenvalues[i];
		for (j=0; j < SVD_MAT_HEIGHT; j++) {
			A[j][i] *= w;
		}
	}

	// Step 2: inverse(A) = V * transpose(A)
	for (i = 0; i < SVD_MAT_WIDTH; i++) {			// rows of Eigenbasis
		for (j = 0; j < SVD_MAT_HEIGHT; j++) {		// colums of A transpose
			InverseA[i][j] = 0.0f;
			for (k = 0; k < SVD_MAT_WIDTH; k++) {	// dimension of overlap
				// A is stored as a transpose so we must flip indices accordingly
				// In other words A[row][column] becomes A[column][row]
				InverseA[i][j] += Eigenbasis[i][k] * A[j][k];
			}
		}
	}

	// compute optimized A and B values
	float frA[4], fgA[4], fbA[4], frB[4], fgB[4], fbB[4];
	int outRed, outGreen, outBlue, repIdx;
	for (j = 0; j < SVD_MAT_WIDTH; j+=2) {
		repIdx = (j>>1);	// repIdx = j / 2
		frA[repIdx] = fgA[repIdx] = fbA[repIdx] = frB[repIdx] = fgB[repIdx] = fbB[repIdx] = 0.0f;
		for (i = 0; i < SVD_MAT_HEIGHT; i++) {
			frA[repIdx] += InverseA[j][i] * red[i];
			fgA[repIdx] += InverseA[j][i] * green[i];
			fbA[repIdx] += InverseA[j][i] * blue[i];
			frB[repIdx] += InverseA[j+1][i] * red[i];
			fgB[repIdx] += InverseA[j+1][i] * green[i];
			fbB[repIdx] += InverseA[j+1][i] * blue[i];
		}
	}

#ifdef USE_PIXEL_UPDATE
	// Calculate new representative colors by using a "fix the error" approach.
	// We multiply the pseudo-inverse weight matrix by the difference between the
	// current and original images to get an "update" that is applied to each
	// representative value. We also limit how big the update can be to avoid
	// colors flying out of bounds.
	unsigned int oldColor;
	// Do top left for 2x2 case or only one for 1x1 case
	tIdx = (width * thready) + threadx;
	oldColor = oldA[tIdx];
	outRed = CLAMP((int)frA[0], -16, 16);
	outGreen = CLAMP((int)fgA[0], -16, 16);
	outBlue = CLAMP((int)fbA[0], -16, 16);
	candidateA[tIdx] = MAKE_ARGB(
		CLAMP((int)MAKE_RED_PIXEL(oldColor) + outRed, 0, 255),
		CLAMP((int)MAKE_GREEN_PIXEL(oldColor) + outGreen, 0, 255),
		CLAMP((int)MAKE_BLUE_PIXEL(oldColor) + outBlue, 0, 255));
	oldA[tIdx] = candidateA[tIdx];

	oldColor = oldB[tIdx];
	outRed = CLAMP((int)frB[0], -16, 16);
	outGreen = CLAMP((int)fgB[0], -16, 16);
	outBlue = CLAMP((int)fbB[0], -16, 16);
	candidateB[tIdx] = MAKE_ARGB(
		CLAMP((int)MAKE_RED_PIXEL(oldColor) + outRed, 0, 255),
		CLAMP((int)MAKE_GREEN_PIXEL(oldColor) + outGreen, 0, 255),
		CLAMP((int)MAKE_BLUE_PIXEL(oldColor) + outBlue, 0, 255));
	oldB[tIdx] = candidateB[tIdx];
#ifdef TWO_BY_TWO
	// Top right case for 2x2
	tIdx = (width * thready) + threadx + 1;
	oldColor = oldA[tIdx];
	outRed = CLAMP((int)frA[1], -16, 16);
	outGreen = CLAMP((int)fgA[1], -16, 16);
	outBlue = CLAMP((int)fbA[1], -16, 16);
	candidateA[tIdx] = MAKE_ARGB(
		CLAMP((int)MAKE_RED_PIXEL(oldColor) + outRed, 0, 255),
		CLAMP((int)MAKE_GREEN_PIXEL(oldColor) + outGreen, 0, 255),
		CLAMP((int)MAKE_BLUE_PIXEL(oldColor) + outBlue, 0, 255));
	oldA[tIdx] = candidateA[tIdx];

	oldColor = oldB[tIdx];
	outRed = CLAMP((int)frB[1], -16, 16);
	outGreen = CLAMP((int)fgB[1], -16, 16);
	outBlue = CLAMP((int)fbB[1], -16, 16);
	candidateB[tIdx] = MAKE_ARGB(
		CLAMP((int)MAKE_RED_PIXEL(oldColor) + outRed, 0, 255),
		CLAMP((int)MAKE_GREEN_PIXEL(oldColor) + outGreen, 0, 255),
		CLAMP((int)MAKE_BLUE_PIXEL(oldColor) + outBlue, 0, 255));
	oldB[tIdx] = candidateB[tIdx];

	// Bottom left case
	tIdx = (width * (thready+1)) + threadx;
	oldColor = oldA[tIdx];
	outRed = CLAMP((int)frA[2], -16, 16);
	outGreen = CLAMP((int)fgA[2], -16, 16);
	outBlue = CLAMP((int)fbA[2], -16, 16);
	candidateA[tIdx] = MAKE_ARGB(
		CLAMP((int)MAKE_RED_PIXEL(oldColor) + outRed, 0, 255),
		CLAMP((int)MAKE_GREEN_PIXEL(oldColor) + outGreen, 0, 255),
		CLAMP((int)MAKE_BLUE_PIXEL(oldColor) + outBlue, 0, 255));
	oldA[tIdx] = candidateA[tIdx];

	oldColor = oldB[tIdx];
	outRed = CLAMP((int)frB[2], -16, 16);
	outGreen = CLAMP((int)fgB[2], -16, 16);
	outBlue = CLAMP((int)fbB[2], -16, 16);
	candidateB[tIdx] = MAKE_ARGB(
		CLAMP((int)MAKE_RED_PIXEL(oldColor) + outRed, 0, 255),
		CLAMP((int)MAKE_GREEN_PIXEL(oldColor) + outGreen, 0, 255),
		CLAMP((int)MAKE_BLUE_PIXEL(oldColor) + outBlue, 0, 255));
	oldB[tIdx] = candidateB[tIdx];

	// Bottom right case
	tIdx = (width * (thready+1)) + threadx + 1;
	oldColor = oldA[tIdx];
	outRed = CLAMP((int)frA[3], -16, 16);
	outGreen = CLAMP((int)fgA[3], -16, 16);
	outBlue = CLAMP((int)fbA[3], -16, 16);
	candidateA[tIdx] = MAKE_ARGB(
		CLAMP((int)MAKE_RED_PIXEL(oldColor) + outRed, 0, 255),
		CLAMP((int)MAKE_GREEN_PIXEL(oldColor) + outGreen, 0, 255),
		CLAMP((int)MAKE_BLUE_PIXEL(oldColor) + outBlue, 0, 255));
	oldA[tIdx] = candidateA[tIdx];

	oldColor = oldB[tIdx];
	outRed = CLAMP((int)frB[3], -16, 16);
	outGreen = CLAMP((int)fgB[3], -16, 16);
	outBlue = CLAMP((int)fbB[3], -16, 16);
	candidateB[tIdx] = MAKE_ARGB(
		CLAMP((int)MAKE_RED_PIXEL(oldColor) + outRed, 0, 255),
		CLAMP((int)MAKE_GREEN_PIXEL(oldColor) + outGreen, 0, 255),
		CLAMP((int)MAKE_BLUE_PIXEL(oldColor) + outBlue, 0, 255));
	oldB[tIdx] = candidateB[tIdx];
#endif // TWO_BY_TWO
#else
	// WARNING - THIS DOES NOT WORK
	// Calculate new representative colors by computing completely new colors.
	// We multiply the pseudo-inverse weight matrix by the actual colors of 
	// the original pixels in our window to get new representative colors.
	// do A candidate. If the new computed color is out of bounds, we discard
	// it and use the old one.
	if ((frA > 255) || (fgA > 255) || (fbA > 255) || (frA < 0) || (fgA < 0) || (fbA < 0)) {
		outRed = MAKE_RED_PIXEL(oldA[tIdx]);
		outGreen = MAKE_GREEN_PIXEL(oldA[tIdx]);
		outBlue = MAKE_BLUE_PIXEL(oldA[tIdx]);
	} else {
		outRed = CLAMP(__float2uint_rn(frA), 0, 255);
		outGreen = CLAMP(__float2uint_rn(fgA), 0, 255);
		outBlue = CLAMP(__float2uint_rn(fbA), 0, 255);
	}
	candidateA[tIdx] = MAKE_ARGB(outRed, outGreen, outBlue);
	oldA[tIdx] = MAKE_ARGB(outRed, outGreen, outBlue);

	// do B candidate
	if ((frB > 255) || (fgB > 255) || (fbB > 255) || (frB < 0) || (fgB < 0) || (fbB < 0)) {
		outRed = MAKE_RED_PIXEL(oldB[tIdx]);
		outGreen = MAKE_GREEN_PIXEL(oldB[tIdx]);
		outBlue = MAKE_BLUE_PIXEL(oldB[tIdx]);
	} else {
		outRed = CLAMP(__float2uint_rn(frB), 0, 255);
		outGreen = CLAMP(__float2uint_rn(fgB), 0, 255);
		outBlue = CLAMP(__float2uint_rn(fbB), 0, 255);
	}
	candidateB[tIdx] = MAKE_ARGB(outRed, outGreen, outBlue);
	oldB[tIdx] = MAKE_ARGB(outRed, outGreen, outBlue);
#endif // USE_PIXEL_UPDATE
}
//////////////////////////////////////////////////////////////////////////////////
//								  Optimization									//
//////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////////
//								 Modulation Bits								//
//////////////////////////////////////////////////////////////////////////////////
// compute modulation bits based on reconstruction of original image
__global__ void compute_modulation_bits(unsigned int *a, unsigned int *b,
										int *redCurrent, int *greenCurrent, 
										int *blueCurrent, float *mod, int width) {
	unsigned int y = blockDim.y * blockIdx.y + threadIdx.y;
	unsigned int x = blockDim.x * blockIdx.x + threadIdx.x;
	int idx = (width * y) + x;

	unsigned int apixel = a[idx];
	int a_red = MAKE_RED_PIXEL(apixel);
	int a_green = MAKE_GREEN_PIXEL(apixel);
	int a_blue = MAKE_BLUE_PIXEL(apixel);

	unsigned int bpixel = b[idx];
	int b_red = MAKE_RED_PIXEL(bpixel);
	int b_green = MAKE_GREEN_PIXEL(bpixel);
	int b_blue = MAKE_BLUE_PIXEL(bpixel);

	unsigned int opixel = tex2D(origRef, x, y);
	int o_red = MAKE_RED_PIXEL(opixel);
	int o_green = MAKE_GREEN_PIXEL(opixel);
	int o_blue = MAKE_BLUE_PIXEL(opixel);

	int3 delta = make_int3(o_red - a_red, o_green - a_green, o_blue - a_blue);
	int3 diff = make_int3(b_red - a_red, b_green - a_green, b_blue - a_blue);

	float dot = (float)(delta.x*diff.x + delta.y*diff.y + delta.z*diff.z) /
				(float)(diff.x*diff.x + diff.y*diff.y + diff.z*diff.z);
	float modbit;

	if (dot < 0.1875f) {
		modbit = 0.00005f;
	} else if (dot < 0.5f) {
		modbit = THREE_EIGHTHS;
	} else if (dot < 0.8125f) {
		modbit = FIVE_EIGHTHS;
	} else {
		modbit = 0.99995f;
	}
#ifdef USE_PIXEL_UPDATE
	float r_modbit = 1.0f - modbit;
	int c_red = ((float)b_red*modbit + (float)a_red*r_modbit);
	int c_green = ((float)b_green*modbit + (float)a_green*r_modbit);
	int c_blue = ((float)b_blue*modbit + (float)a_blue*r_modbit);

	redCurrent[idx] = o_red - c_red;
	greenCurrent[idx] = o_green - c_green;
	blueCurrent[idx] = o_blue - c_blue;
#else
	// WARNING - THIS DOES NOT WORK
	redCurrent[idx] = o_red;
	greenCurrent[idx] = o_green;
	blueCurrent[idx] = o_blue;
#endif
	mod[idx] = modbit;
}
//////////////////////////////////////////////////////////////////////////////////
//								 Modulation Bits								//
//////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////////
//								A and B prototypes								//
//////////////////////////////////////////////////////////////////////////////////
__global__ void make_a_b_prototypes(unsigned int *a_proto, unsigned int *b_proto,
									int width, int height) {
	unsigned int y = blockDim.y * blockIdx.y + threadIdx.y;
	unsigned int x = blockDim.x * blockIdx.x + threadIdx.x;
	int idx = (width * y) + x;

	// get original, delta, and axis pixel
	int opixel = tex2D(origRef, x, y);
	int fpixel = tex2D(texRef, x, y);

	int orig_red = MAKE_RED_PIXEL(opixel);
	int orig_green = MAKE_GREEN_PIXEL(opixel);
	int orig_blue = MAKE_BLUE_PIXEL(opixel);
	int delta_red = orig_red - MAKE_RED_PIXEL(fpixel);
	int delta_green = orig_green - MAKE_GREEN_PIXEL(fpixel);
	int delta_blue = orig_blue - MAKE_BLUE_PIXEL(fpixel);
	int axis_red = CLAMP(abs(delta_red), 0, 255);
	int axis_green = CLAMP(abs(delta_green), 0, 255);
	int axis_blue = CLAMP(abs(delta_blue), 0, 255);

	// (delta dot axis) / (axis dot axis)
	int dot = ((delta_red * axis_red) + (delta_green * axis_green) + 
			  (delta_blue * axis_blue)) / ((axis_red * axis_red) + 
			  (axis_green * axis_green) + (axis_blue * axis_blue));

	/*if (dot < 0) {
		a_proto[idx] = MAKE_ARGB((CLAMP(orig_red + axis_red * dot, 0, 255)), 
								 (CLAMP(orig_green + axis_green * dot, 0, 255)), 
								 (CLAMP(orig_blue + axis_blue * dot, 0, 255)));

		b_proto[idx] = MAKE_ARGB((CLAMP(orig_red - axis_red * dot, 0, 255)), 
								 (CLAMP(orig_green - axis_green * dot, 0, 255)), 
								 (CLAMP(orig_blue - axis_blue * dot, 0, 255)));
	} else {
		b_proto[idx] = MAKE_ARGB((CLAMP(orig_red + axis_red * dot, 0, 255)), 
								 (CLAMP(orig_green + axis_green * dot, 0, 255)), 
								 (CLAMP(orig_blue + axis_blue * dot, 0, 255)));

		a_proto[idx] = MAKE_ARGB((CLAMP(orig_red - axis_red * dot, 0, 255)), 
								 (CLAMP(orig_green - axis_green * dot, 0, 255)), 
								 (CLAMP(orig_blue - axis_blue * dot, 0, 255)));
	}*/
	b_proto[idx] = MAKE_ARGB((CLAMP(orig_red + axis_red, 0, 255)), 
		(CLAMP(orig_green + axis_green, 0, 255)), 
		(CLAMP(orig_blue + axis_blue, 0, 255)));

	a_proto[idx] = MAKE_ARGB((CLAMP(orig_red - axis_red, 0, 255)), 
		(CLAMP(orig_green - axis_green, 0, 255)), 
		(CLAMP(orig_blue - axis_blue, 0, 255)));
}
//////////////////////////////////////////////////////////////////////////////////
//								A and B prototypes								//
//////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////////
//								bilinear upscaling								//
//////////////////////////////////////////////////////////////////////////////////

// (1-alpha)*A + alpha*B. alpha must be between 0~1.
__device__ float lerp(const float A, const float B, float alpha)
{
	return (((1-alpha)*A) + (alpha*B));
}

__global__ void bilinear_resize4x4(unsigned int *out, int width, int height, bool flag) {
	int x, x1, y, y1;	// x1 and y1 for texture wrapping purposes
	int j = blockDim.x * blockIdx.x + threadIdx.x;
	int i = blockDim.y * blockIdx.y + threadIdx.y;
	unsigned int a, b, c, d;
	float x_diff, y_diff, blue, red, green;

	// do the 4x4 bilinear scaling
	if (flag == true) {
		x = CLAMP((j-2)>>2, 0, width);
		y = CLAMP((i-2)>>2, 0, height);
		x1 = x+1;
		y1 = y+1;
		x_diff = (CLAMP(j-2, 0, width) - (x*4)) * ONE_FOURTH; // j%4 = j - (x*4)
		y_diff = (CLAMP(i-2, 0, height) - (y*4)) * ONE_FOURTH;
	} else {
		x = (j>>2);
		y = (i>>2);
		x1 = CLAMP(x+1, 0, (width>>2)-1);
		y1 = CLAMP(y+1, 0, (width>>2)-1);
		x_diff = (j - (x*4)) * ONE_FOURTH; // j%4 = j - (x*4)
		y_diff = (i - (y*4)) * ONE_FOURTH;
	}

	// fetch the 32bit ARGB unsigned ints for the four pixels
	a = tex2D(texRef, x, y);
	b = tex2D(texRef, x1, y);
	c = tex2D(texRef, x, y1);
	d = tex2D(texRef, x1, y1);

	// blue channel
	blue = lerp(lerp((float)MAKE_BLUE_PIXEL(a), (float)MAKE_BLUE_PIXEL(b), x_diff),
		lerp((float)MAKE_BLUE_PIXEL(c), (float)MAKE_BLUE_PIXEL(d), x_diff),
		y_diff);

	// green channel
	green = lerp(lerp((float)MAKE_GREEN_PIXEL(a), (float)MAKE_GREEN_PIXEL(b), x_diff),
		lerp((float)MAKE_GREEN_PIXEL(c), (float)MAKE_GREEN_PIXEL(d), x_diff),
		y_diff);

	// red channel
	red = lerp(lerp((float)MAKE_RED_PIXEL(a), (float)MAKE_RED_PIXEL(b), x_diff),
		lerp((float)MAKE_RED_PIXEL(c), (float)MAKE_RED_PIXEL(d), x_diff),
		y_diff);

	// store the result
	out[i*width + j] = MAKE_ARGB( CLAMP(__float2uint_rn(red), 0, 255),
		CLAMP(__float2uint_rn(green), 0, 255),
		CLAMP(__float2uint_rn(blue), 0, 255) );
}

__global__ void bilinear_resize8x4(unsigned int *out, int width, int height, bool flag) {
	int x, x1, y, y1;	// x1 and y1 for texture wrapping purposes
	int j = blockDim.x * blockIdx.x + threadIdx.x;
	int i = blockDim.y * blockIdx.y + threadIdx.y;
	unsigned int a, b, c, d;
	float x_diff, y_diff, blue, red, green;

	// do the 8x4 bilinear scaling
	x = j>>3;		// (j / 8)
	y = i>>2;		// (i / 4)
	x1 = (x+1) & ((width>>3) - 1);	// (x+1) % (width / 8)
	y1 = (y+1) & ((height>>2) - 1);	// (y+1) % (height / 4)
	x_diff = (j - (x*8)) * ONE_EIGHTH;	// j%8 = j - (x*8)
	y_diff = (i - (y*4)) * ONE_FOURTH;	// i%4 = j - (y*4)

	// fetch the 32bit ARGB unsigned ints for the four pixels
	a = tex2D(texRef, x, y);
	b = tex2D(texRef, x1, y);
	c = tex2D(texRef, x, y1);
	d = tex2D(texRef, x1, y1);

	// blue channel
	blue = lerp(lerp((float)MAKE_BLUE_PIXEL(a), (float)MAKE_BLUE_PIXEL(b), x_diff),
		lerp((float)MAKE_BLUE_PIXEL(c), (float)MAKE_BLUE_PIXEL(d), x_diff),
		y_diff);

	// green channel
	green = lerp(lerp((float)MAKE_GREEN_PIXEL(a), (float)MAKE_GREEN_PIXEL(b), x_diff),
		lerp((float)MAKE_GREEN_PIXEL(c), (float)MAKE_GREEN_PIXEL(d), x_diff),
		y_diff);

	// red channel
	red = lerp(lerp((float)MAKE_RED_PIXEL(a), (float)MAKE_RED_PIXEL(b), x_diff),
		lerp((float)MAKE_RED_PIXEL(c), (float)MAKE_RED_PIXEL(d), x_diff),
		y_diff);

	// offset by 2x2 to center representative values
	if (flag == true) {
		// since width and height HAVE to be a power of 2, we can take a 
		// shortcut for modulo here. x%y = x & (y-1) when y is a power of 2
		out[((i+2) & (height-1))*width + ((j+2) & (width-1))] = MAKE_ARGB(
			CLAMP(__float2uint_rn(red), 0, 255),
			CLAMP(__float2uint_rn(green), 0, 255),
			CLAMP(__float2uint_rn(blue), 0, 255));
	} else {
		out[i*width + j] = MAKE_ARGB(
			CLAMP(__float2uint_rn(red), 0, 255),
			CLAMP(__float2uint_rn(green), 0, 255),
			CLAMP(__float2uint_rn(blue), 0, 255));
	}
}
//////////////////////////////////////////////////////////////////////////////////
//								 bilinear upscaling								//
//////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////////
//							low-pass wavelet filter								//
//////////////////////////////////////////////////////////////////////////////////
__global__ void linear_wavelet_transform_rows(int width, int height, 
											  unsigned int* out) {
	int n; 
	int y = blockDim.y * blockIdx.y + threadIdx.y;
	int k = blockDim.x * blockIdx.x + threadIdx.x;
	unsigned int pixel;
	float filter;

	float3 s = make_float3(0.0f, 0.0f, 0.0f);

	// we start out negative because we want to center the filter on
	// the current pixel
	for (int m = 0 - HALF_FILTER_LENGTH; m < HALF_FILTER_LENGTH; m++) {
		n = 2 * k + m;
		if (n < 0) {
			n = 0 - n;
		}
		if (n >= width) {
			n -= 2 * (1 + n - width);
		}
		pixel = tex2D(texRef, CLAMP(n, 0, width-1), CLAMP(y, 0, height-1));
		filter = tex1Dfetch(filterRef, (int)(m + HALF_FILTER_LENGTH));

		s.x += filter * __uint2float_rn(MAKE_RED_PIXEL(pixel));
		s.y += filter * __uint2float_rn(MAKE_GREEN_PIXEL(pixel));
		s.z += filter * __uint2float_rn(MAKE_BLUE_PIXEL(pixel));
	}

	out[y*width + k] = MAKE_ARGB(CLAMP(__float2uint_rn(s.x), 0, 255),
								 CLAMP(__float2uint_rn(s.y), 0, 255),
								 CLAMP(__float2uint_rn(s.z), 0, 255));
}

__global__ void linear_wavelet_transform_cols(int width, int height, 
											  unsigned int* out) {
	int n; 
	int k = blockDim.y * blockIdx.y + threadIdx.y;
	int x = blockDim.x * blockIdx.x + threadIdx.x;
	unsigned int pixel;
	float filter;

	float3 s = make_float3(0.0f, 0.0f, 0.0f);

	// we start out negative because we want to center the filter on
	// the current pixel
	for (int m = 0 - HALF_FILTER_LENGTH; m < HALF_FILTER_LENGTH; m++) {
		n = 2 * k + m;
		if (n < 0) {
			n = 0 - n;
		}
		if (n >= height) {
			n -= 2 * (1 + n - height);
		}
		pixel = tex2D(texRef, CLAMP(x, 0, width-1), CLAMP(n, 0, height-1));
		filter = tex1Dfetch(filterRef, (int)(m + HALF_FILTER_LENGTH));

		s.x += filter * __uint2float_rn(MAKE_RED_PIXEL(pixel));
		s.y += filter * __uint2float_rn(MAKE_GREEN_PIXEL(pixel));
		s.z += filter * __uint2float_rn(MAKE_BLUE_PIXEL(pixel));
	}

	out[k*width + x] = MAKE_ARGB(CLAMP(__float2uint_rn(s.x), 0, 255),
								 CLAMP(__float2uint_rn(s.y), 0, 255),
								 CLAMP(__float2uint_rn(s.z), 0, 255));
}

//////////////////////////////////////////////////////////////////////////////////
//	Calculates the low-pass wavelet transform of the input data and stores the	//
//	result in the same location. First transforms the rows and then the colums	//
//	of the input image.															//
//																				//
//	wavelet		Initial input data. The result is also stored here				//
//	temp		Pointer to scratch space in memory that is the same size as the	//
//				input data														//
//	num			The number of times the low-pass filter is applied				//
//																				//
//////////////////////////////////////////////////////////////////////////////////
void linear_wavelet_transform(unsigned int * wavelet, unsigned int *temp, int num,
							  textureReference *texRefPtr, int width, int height, 
							  int scan_width, cudaChannelFormatDesc *channelDesc,
							  dim3 default_block){
	size_t offset;
	int i;
	for (i = 1; i <= num; i++) {
		// bind the input data to the texture
		cutilSafeCall(cudaBindTexture2D(&offset, texRefPtr, (const void*)wavelet, 
			channelDesc, width, height, scan_width));

		// calulate the grid dimensions and call the kernel
		dim3 grid0(width / (BLOCK_WIDTH * 2 * i), height / (BLOCK_HEIGHT * i));
		linear_wavelet_transform_rows<<<grid0, default_block>>> (width, height, temp);

		// bind the new data to the texture
		cutilSafeCall(cudaBindTexture2D(&offset, texRefPtr, (const void*)temp, 
			channelDesc, width, height, scan_width));
		
		// calulate the grid dimensions and call the kernel
		dim3 grid1(width / (BLOCK_WIDTH * i), height / (BLOCK_HEIGHT * 2 * i));
		linear_wavelet_transform_cols<<<grid1, default_block>>> (width, height, wavelet);
	}
}
//////////////////////////////////////////////////////////////////////////////////
//							low-pass wavelet filter								//
//////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////////
//	External method to call CUDA kernels for compression. Places the compressed	//
//	data in	h_out.																//
//																				//
//	width		Width of the input image										//
//	height		Height of the input image										//
//	scan_width	Size in bytes of a single row of pixels							//
//	h_in		Pointer to the input buffer on the host							//
//	h_out		Pointer to the output buffer on the	host						//
//																				//
//////////////////////////////////////////////////////////////////////////////////
extern "C" int cuda_pvr_compress(int width, int height, int scan_width, 
								  unsigned char* h_in, unsigned char* h_out,
								  int num_optimizations) {
	// image dimensions must be powers of 2
	if ((width & (width - 1)) != 0 || (height & (height - 1)) != 0) {
		printf("Error: Image dimensions must be powers of 2.\n");
		printf("Aborting compression...\n");

		return 1;
	}

	// declare device memory pointers and block dimensions
	unsigned int *d_bits,*d_out,*d_temp,*d_wavelet,*d_axis,*d_aproto,*d_bproto;
	int *d_redCurrent, *d_greenCurrent, *d_blueCurrent, *d_err;
	float *d_filter, *d_mod;
#ifdef TWO_BY_TWO
    float *d_svdMatrixTL, *d_svdMatrixTR, *d_svdMatrixBL, *d_svdMatrixBR;
#else
    float *d_svdMatrix;
#endif // TWO_BY_TWO
	int i, tex_size = height * scan_width;
	dim3 default_block(BLOCK_WIDTH, BLOCK_HEIGHT);
	dim3 default_grid(width / (BLOCK_WIDTH), height / (BLOCK_HEIGHT));
	dim3 quarterGrid(default_grid.x>>2, default_grid.y>>2);
    dim3 svdGrid(default_grid.x * SVD_FACTOR, default_grid.y * SVD_FACTOR);
	size_t offset;

	//////////////////////////////////////////////////////////////////////////////
	//////////////////////////////////////////////////////////////////////////////

	// initialize cuda and allocate device memory
	CUT_DEVICE_INIT(1, "");
	cutilSafeCall(cudaMalloc((void **)&d_bits, tex_size));
	cutilSafeCall(cudaMalloc((void **)&d_temp, tex_size));
	cutilSafeCall(cudaMalloc((void **)&d_wavelet, tex_size));
	cutilSafeCall(cudaMalloc((void **)&d_axis, tex_size));
	cutilSafeCall(cudaMalloc((void **)&d_aproto, tex_size));
	cutilSafeCall(cudaMalloc((void **)&d_bproto, tex_size));
	cutilSafeCall(cudaMalloc((void **)&d_redCurrent, width * height * sizeof(int)));
	cutilSafeCall(cudaMalloc((void **)&d_greenCurrent, width * height * sizeof(int)));
	cutilSafeCall(cudaMalloc((void **)&d_blueCurrent, width * height * sizeof(int)));
	cutilSafeCall(cudaMalloc((void **)&d_err, sizeof(int)));
	cutilSafeCall(cudaMalloc((void **)&d_filter, FILTER_LENGTH * sizeof(float)));
	cutilSafeCall(cudaMalloc((void **)&d_mod, width * height * sizeof(float)));
	cutilSafeCall(cudaMalloc((void **)&d_out, 2 * (width>>2) * (height>>2) * 
		sizeof(unsigned int)));

	cutilSafeCall(cudaMemcpy(d_bits, h_in, tex_size, cudaMemcpyHostToDevice));
	cutilSafeCall(cudaMemcpy(d_wavelet, d_bits, tex_size, cudaMemcpyDeviceToDevice));
	cutilSafeCall(cudaMemcpy(d_filter, wavelet_filter, FILTER_LENGTH * sizeof(float), 
							 cudaMemcpyHostToDevice));
#ifdef TWO_BY_TWO
    cutilSafeCall(cudaMalloc((void **)&d_svdMatrixTL, SVD_MAT_HEIGHT * sizeof(float)));
    cutilSafeCall(cudaMalloc((void **)&d_svdMatrixTR, SVD_MAT_HEIGHT * sizeof(float)));
    cutilSafeCall(cudaMalloc((void **)&d_svdMatrixBL, SVD_MAT_HEIGHT * sizeof(float)));
    cutilSafeCall(cudaMalloc((void **)&d_svdMatrixBR, SVD_MAT_HEIGHT * sizeof(float)));
    cutilSafeCall(cudaMemcpy(d_svdMatrixTL, MwTL, SVD_MAT_HEIGHT * sizeof(float), 
							 cudaMemcpyHostToDevice));
    cutilSafeCall(cudaMemcpy(d_svdMatrixTR, MwTR, SVD_MAT_HEIGHT * sizeof(float), 
							 cudaMemcpyHostToDevice));
    cutilSafeCall(cudaMemcpy(d_svdMatrixBL, MwBL, SVD_MAT_HEIGHT * sizeof(float), 
							 cudaMemcpyHostToDevice));
    cutilSafeCall(cudaMemcpy(d_svdMatrixBR, MwBR, SVD_MAT_HEIGHT * sizeof(float), 
							 cudaMemcpyHostToDevice));
#else
    cutilSafeCall(cudaMalloc((void **)&d_svdMatrix, SVD_MAT_HEIGHT * sizeof(float)));
	cutilSafeCall(cudaMemcpy(d_svdMatrix, Mw, SVD_MAT_HEIGHT * sizeof(float), 
							 cudaMemcpyHostToDevice));
#endif // TWO_BY_TWO

	//////////////////////////////////////////////////////////////////////////////
	//////////////////////////////////////////////////////////////////////////////

	// map the input buffer as texture memory
	cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<unsigned int>();
	cudaChannelFormatDesc channelDesc2 = cudaCreateChannelDesc<float>();
	cudaChannelFormatDesc channelDesc3 = cudaCreateChannelDesc<int>();

	textureReference *texRefPtr,*origRefPtr,*filterRefPtr,*modRefPtr,
		*redRefPtr,*greenRefPtr,*blueRefPtr; 

	cudaGetTextureReference((const textureReference**)&texRefPtr, "texRef");
	cudaGetTextureReference((const textureReference**)&origRefPtr, "origRef");
	cudaGetTextureReference((const textureReference**)&filterRefPtr, "filterRef");
	cudaGetTextureReference((const textureReference**)&modRefPtr, "modRef");
	cudaGetTextureReference((const textureReference**)&redRefPtr, "redCurrentRef");
	cudaGetTextureReference((const textureReference**)&greenRefPtr, "greenCurrentRef");
	cudaGetTextureReference((const textureReference**)&blueRefPtr, "blueCurrentRef");
    
#ifdef TWO_BY_TWO
    textureReference *svdRefPtrTL, *svdRefPtrTR, *svdRefPtrBL, *svdRefPtrBR;
    cudaGetTextureReference((const textureReference**)&svdRefPtrTL, "svdMatRefTL");
    cudaGetTextureReference((const textureReference**)&svdRefPtrTR, "svdMatRefTR");
    cudaGetTextureReference((const textureReference**)&svdRefPtrBL, "svdMatRefBL");
    cudaGetTextureReference((const textureReference**)&svdRefPtrBR, "svdMatRefBR");
#else
    textureReference *svdRefPtr;
    cudaGetTextureReference((const textureReference**)&svdRefPtr, "svdMatRef");
#endif // TWO_BY_TWO

	// set the addressing and filter modes for all 2D textures
	texRefPtr->addressMode[0] = cudaAddressModeClamp;
	texRefPtr->addressMode[1] = cudaAddressModeClamp;
	texRefPtr->addressMode[2] = cudaAddressModeClamp;
	texRefPtr->filterMode = cudaFilterModePoint;

	origRefPtr->addressMode[0] = cudaAddressModeClamp;
	origRefPtr->addressMode[1] = cudaAddressModeClamp;
	origRefPtr->addressMode[2] = cudaAddressModeClamp;
	origRefPtr->filterMode = cudaFilterModePoint;

	modRefPtr->addressMode[0] = cudaAddressModeClamp;
	modRefPtr->addressMode[1] = cudaAddressModeClamp;
	modRefPtr->addressMode[2] = cudaAddressModeClamp;
	modRefPtr->filterMode = cudaFilterModePoint;

	redRefPtr->addressMode[0] = cudaAddressModeClamp;
	redRefPtr->addressMode[1] = cudaAddressModeClamp;
	redRefPtr->addressMode[2] = cudaAddressModeClamp;
	redRefPtr->filterMode = cudaFilterModePoint;

	blueRefPtr->addressMode[0] = cudaAddressModeClamp;
	blueRefPtr->addressMode[1] = cudaAddressModeClamp;
	blueRefPtr->addressMode[2] = cudaAddressModeClamp;
	blueRefPtr->filterMode = cudaFilterModePoint;

	greenRefPtr->addressMode[0] = cudaAddressModeClamp;
	greenRefPtr->addressMode[1] = cudaAddressModeClamp;
	greenRefPtr->addressMode[2] = cudaAddressModeClamp;
	greenRefPtr->filterMode = cudaFilterModePoint;

	// bind textures
	cutilSafeCall(cudaBindTexture2D(&offset, origRefPtr, (const void*)d_bits, 
		&channelDesc, width, height, scan_width));
	cutilSafeCall(cudaBindTexture(&offset, filterRefPtr, (const void*)d_filter, 
		&channelDesc2, FILTER_LENGTH * sizeof(float)));
	cutilSafeCall(cudaBindTexture2D(&offset, modRefPtr, (const void*)d_mod, 
		&channelDesc2, width, height, width * sizeof(float)));
	cutilSafeCall(cudaBindTexture2D(&offset, redRefPtr, (const void*)d_redCurrent, 
		&channelDesc3, width, height, width * sizeof(int)));
	cutilSafeCall(cudaBindTexture2D(&offset, greenRefPtr, (const void*)d_greenCurrent, 
		&channelDesc3, width, height, width * sizeof(int)));
	cutilSafeCall(cudaBindTexture2D(&offset, blueRefPtr, (const void*)d_blueCurrent, 
		&channelDesc3, width, height, width * sizeof(int)));
#ifdef TWO_BY_TWO
    cutilSafeCall(cudaBindTexture(&offset, svdRefPtrTL, (const void*)d_svdMatrixTL, 
                                  &channelDesc2, SVD_MAT_HEIGHT * sizeof(float)));
    cutilSafeCall(cudaBindTexture(&offset, svdRefPtrTR, (const void*)d_svdMatrixTR, 
                                  &channelDesc2, SVD_MAT_HEIGHT * sizeof(float)));
    cutilSafeCall(cudaBindTexture(&offset, svdRefPtrBL, (const void*)d_svdMatrixBL, 
                                  &channelDesc2, SVD_MAT_HEIGHT * sizeof(float)));
    cutilSafeCall(cudaBindTexture(&offset, svdRefPtrBR, (const void*)d_svdMatrixBR, 
                                  &channelDesc2, SVD_MAT_HEIGHT * sizeof(float)));
#else
    cutilSafeCall(cudaBindTexture(&offset, svdRefPtr, (const void*)d_svdMatrix, 
                                  &channelDesc2, SVD_MAT_HEIGHT * sizeof(float)));
#endif // TWO_BY_TWO

	//////////////////////////////////////////////////////////////////////////////
	//////////////////////////////////////////////////////////////////////////////

	// apply low-pass filter to initial image
	linear_wavelet_transform(d_wavelet, d_temp, 2, texRefPtr, width, height, 
							 scan_width, &channelDesc, default_block);
	
	// 4x4 bilinear upscale
	cutilSafeCall(cudaBindTexture2D(&offset, texRefPtr, (const void*)d_wavelet, 
				  &channelDesc, width, height, scan_width));
	bilinear_resize4x4<<<default_grid, default_block>>> (d_temp, width, height, false);
	swap((void**)&d_temp, (void**)&d_wavelet);

	//////////////////////////////////////////////////////////////////////////////
	//////////////////////////////////////////////////////////////////////////////

	// make the inital A and B prototypes
	cutilSafeCall(cudaBindTexture2D(&offset, texRefPtr, (const void*)d_wavelet, 
				  &channelDesc, width, height, scan_width));
	make_a_b_prototypes<<<default_grid, default_block>>> (d_aproto, d_bproto, 
		width, height);

	// apply low-pass filter to A and B prototypes to get initial candidates
	linear_wavelet_transform(d_aproto, d_temp, 2, texRefPtr, width, height, 
		scan_width, &channelDesc, default_block);
	linear_wavelet_transform(d_bproto, d_temp, 2, texRefPtr, width, height, 
		scan_width, &channelDesc, default_block);

	// copy the prototypes into temporary storage (re-using d_axis and d_wavelet
	// for temporary storage)
	cudaMemcpy(d_axis, d_aproto, tex_size, cudaMemcpyDeviceToDevice);
	cudaMemcpy(d_wavelet, d_bproto, tex_size, cudaMemcpyDeviceToDevice);

	//////////////////////////////////////////////////////////////////////////////
	//////////////////////////////////////////////////////////////////////////////

	// optimize A and B candidates and encode final image
	for (i = 0; i < num_optimizations; i++) {
		// expand the A and B candidates to full resolution
		cutilSafeCall(cudaBindTexture2D(&offset, texRefPtr, (const void*)d_aproto, 
			&channelDesc, width, height, scan_width));
		bilinear_resize4x4<<<default_grid, default_block>>> (d_temp, width, height, true);
		swap((void**)&d_temp, (void**)&d_aproto);

		cutilSafeCall(cudaBindTexture2D(&offset, texRefPtr, (const void*)d_bproto, 
			&channelDesc, width, height, scan_width));
		bilinear_resize4x4<<<default_grid, default_block>>> (d_temp, width, height, true);
		swap((void**)&d_temp, (void**)&d_bproto);

		// get the modulation bits 
		compute_modulation_bits<<<default_grid, default_block>>> (d_aproto, d_bproto, 
			d_redCurrent, d_greenCurrent, d_blueCurrent, d_mod, width);

		// make optimization call
		cutilSafeCall(cudaBindTexture2D(&offset, texRefPtr, (const void*)d_bits, 
										&channelDesc, width, height, scan_width));
#ifdef USE_SVD
		svd_optimize<<<svdGrid, default_block>>>(d_aproto, d_bproto, d_axis,
			d_wavelet, width, height, d_err);
#else
		moore_penrose_optimize<<<quarterGrid, default_block>>>(d_aproto, d_bproto, d_axis,
			d_wavelet, width, height, d_err);
#endif // USE_SVD
	}

	// check for SVD error
	int kernelError;
	cudaMemcpy(&kernelError, d_err, sizeof(int), cudaMemcpyDeviceToHost);
	if (kernelError == 1) {
		printf("One or more SVD matrices did not converge in 30 iterations.\n");
		printf("Compression failed.\n");
		printf("Please try again.\n\n");
		printf("Press ENTER to exit.\n");
		getchar();
		exit(EXIT_FAILURE);
	}

	//////////////////////////////////////////////////////////////////////////////
	//////////////////////////////////////////////////////////////////////////////

#ifdef DECOMPRESS_PVR
	// expand the A and B candidates to full resolution
	cutilSafeCall(cudaBindTexture2D(&offset, texRefPtr, (const void*)d_aproto, 
		&channelDesc, width, height, scan_width));
	bilinear_resize4x4<<<default_grid, default_block>>> (d_temp, width, height, true);
	swap((void**)&d_temp, (void**)&d_aproto);

	cutilSafeCall(cudaBindTexture2D(&offset, texRefPtr, (const void*)d_bproto, 
		&channelDesc, width, height, scan_width));
	bilinear_resize4x4<<<default_grid, default_block>>> (d_temp, width, height, true);
	swap((void**)&d_temp, (void**)&d_bproto);

	decompress<<<default_grid, default_block>>>(d_aproto, d_bproto, d_wavelet, d_mod, width);
	cutilSafeCall(cudaMemcpy(h_out, d_wavelet, tex_size, cudaMemcpyDeviceToHost));
#else

	// return compressed data
	encode_texture<<<quarterGrid, default_block>>>(d_aproto, d_bproto, d_mod, width, 
		false, d_out, width>>2, height>>2);
	cutilSafeCall(cudaMemcpy(h_out, d_out, 2 * (width>>2) * (height>>2) * 
		sizeof(unsigned int), cudaMemcpyDeviceToHost));
	//cutilSafeCall(cudaMemcpy(h_out, d_aproto, tex_size, cudaMemcpyDeviceToHost));

#ifdef GET_RMS_ERROR
	// expand the A and B candidates to full resolution
	cutilSafeCall(cudaBindTexture2D(&offset, texRefPtr, (const void*)d_aproto, 
		&channelDesc, width, height, scan_width));
	bilinear_resize4x4<<<default_grid, default_block>>> (d_temp, width, height, true);
	swap((void**)&d_temp, (void**)&d_aproto);

	cutilSafeCall(cudaBindTexture2D(&offset, texRefPtr, (const void*)d_bproto, 
		&channelDesc, width, height, scan_width));
	bilinear_resize4x4<<<default_grid, default_block>>> (d_temp, width, height, true);
	swap((void**)&d_temp, (void**)&d_bproto);

	// get the per pixel error
	unsigned int *h_rms, *d_rms, j;
	float error, mean_error, total_error = 0.0f;
	h_rms = (unsigned int*)malloc(width * height * sizeof(int));
	cutilSafeCall(cudaMalloc((void **)&d_rms, width * height * sizeof(int)));

	rms_error<<<default_grid, default_block>>>(d_bits, d_aproto, d_bproto, d_rms, width);
	cutilSafeCall(cudaMemcpy(h_rms, d_rms, width * height * sizeof(int), 
		cudaMemcpyDeviceToHost));
	cudaThreadSynchronize();
	for (i = 0; i < height; i++) {
		for (j = 0; j < width; j++) {
			total_error += h_rms[i*width + j];
		}
	}
	mean_error = total_error / (width * height * 3);
	error = sqrtf(mean_error);
	printf("RMS Error is %.3f.\n", error);
#endif // GET_RMS_ERROR
#endif // DECOMPRESS_PVR

	// cleanup
	cudaThreadExit();

	return 0;
}